Tree topology diagnostics

Long-billed hermit song cultural evolution

Source code and data found at https://github.com/maRce10/lbh_cultural_evolution_revBayes

Load packages

Custom functions and global options

knitr::opts_chunk$set(out.width = "100%", out.height = "100%", echo = TRUE) 

source('~/Dropbox/Projects/lbh_cultural_evolution_revBayes/source/custom_treespace.R')

theme_set(theme_classic(base_size = 10))

fill_color <- viridis::mako(10, alpha = 0.5)[7]

iterations <- 5000
chains <- 4
priors <- c(set_prior("normal(0, 1.5)", class = "Intercept"), set_prior("normal(0, 1.5)", class = "b"))

# to create several posterior predictive check plots out of a brms fit
custom_ppc1 <- function(fit, group = NULL, ndraws = 30) {
  
    ppc_dens <-
        pp_check(fit,
                 ndraws = ndraws,
                 type = 'dens_overlay_grouped',
                 group = group)  # shows dens_overlay plot by default
    pp_mean <-
        pp_check(
            fit,
            type = "stat_grouped",
            stat = "mean",
            group = group,
            ndraws = ndraws
        )
    pp_scat <-
        pp_check(fit,
                 type = "scatter_avg_grouped",
                 group = group,
                 ndraws = ndraws)
    pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
    pp_loo <- pp_check(fit, type = "loo_pit_qq", ndraws = ndraws)
    pp_error <-
        pp_check(fit,
                 type = "error_scatter_avg_grouped",
                 ndraws = ndraws,
                 group = group)
    plot_l <-
        list(ppc_dens, pp_mean, pp_scat, pp_error,  pp_stat2, pp_loo)
    
    plot_l[c(1, 3:length(plot_l))] <-
        lapply(plot_l[c(1, 3:length(plot_l))], function(x)
            x  + scale_color_viridis_d(
                begin = 0.1,
                end = 0.8,
                alpha = 0.5
            ) + scale_fill_viridis_d(
                begin = 0.1,
                end = 0.8,
                alpha = 0.5
            ) + theme_classic())     

    print(plot_grid(plotlist = plot_l, ncol = 2))
}

custom_ppc <- function(...) suppressMessages(custom_ppc1(...))

Read data

rb.trees <- readRDS("./output/revbayes_output_in_single_R_object.RDS")

rb.trees <- rb.trees[grep("run_", names(rb.trees), invert = TRUE, value = TRUE)]

trees_diags <- data.frame(model = names(rb.trees))

# get lek name
trees_diags$lek <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 1)

# substitution model
trees_diags$subs <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 2)

# period
trees_diags$period <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 3)

# and which fossils were used
trees_diags$fossils <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 4)

trees_diags$align <- sapply(strsplit(trees_diags$model, "_", fixed = TRUE), "[", 5)

trees_diags$n_trees <- sapply(rb.trees, function(x) length(x$trees))

# Number of models by iterations and lek
kbl <-
    knitr::kable(as.matrix(table(trees_diags$lek, trees_diags$n_trees)), caption = "Number of models with a specific number of trees by lek")

kableExtra::kable_styling(kbl)

# Number of models by iterations and lek
kbl <-
    knitr::kable(as.matrix(table(trees_diags$lek, trees_diags$subs)), caption = "Number of models with a specific number of trees by lek")

kableExtra::kable_styling(kbl)

1 Topological space estimation

• Pairwise topological distance between all trees for each lek
• 100 trees for each model evenly spaced along the chain
• Comparing based on shared song types (pairwise shared tips)
• Some trees were removed if NAs were produced when comparing topologies (i.e. non-comparable topologies)

1.1 All fossils

#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")

tree_names <- grep("early", names(rb.trees), invert = TRUE, value = TRUE)

pnts_lks <- lapply(sel_leks, function(i) {
  
  print(i)
  lek_trees <- grep(i, tree_names, value = TRUE)
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))
  
  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "RF", cl = 10), silent = TRUE)

  if (!is(pnts, "try-error")) 
    return(pnts) else
      return(NULL)
})

names(pnts_lks) <- sel_leks

sapply(pnts_lks, is.null)

saveRDS(pnts_lks, file = "./output/topology/tree_distance_all_fossils_pairwise_shared.RDS")

1.2 Early fossils

#selected leks
sel_leks <- c("SUR", "CCE", "HC1", "BR1", "TR1")
# sel_leks <- c("SUR", "TR1")

tree_names <- grep("early", names(rb.trees), value = TRUE)

pnts_lks <- lapply(sel_leks, function(i) {
  
  print(i)
  lek_trees <- grep(i, tree_names, value = TRUE)
  
  multi_list <- lapply(lek_trees, function(x) read.tree(file.path("./output/most_recent_revbayes_models/", x), keep.multi = TRUE))
  
  names(multi_list) <- sapply(basename(lek_trees), function(x) paste(strsplit(x, "_", fixed = TRUE)[[1]][1:5], collapse = "_"))
  
  pnts <- try(custom_treespace(multi_list, n.points = 100, method = "RF", cl = 10), silent = TRUE)

  if (!is(pnts, "try-error")) 
    return(pnts) else
      return(NULL)
})

names(pnts_lks) <- sel_leks

sapply(pnts_lks, is.null)

saveRDS(pnts_lks, file = "./output/topology/tree_distance_early_fossils_pairwise_shared.RDS")

2 Plot topological spaces

2.1 All fossils

2.1.1 New and old data sets

pnts_lks <- readRDS("./output/topology/tree_distance_all_fossils_pairwise_shared.RDS")

ggs <- lapply(pnts_lks, function(X) {

    ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, ncol = 2) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") + 
      ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain)))) 
    })

ggs

$SUR

$CCE

$HC1

$BR1

$TR1

2.1.2 Old data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  
  if (nrow(X) > 0)
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))
    })

ggs

$SUR

$CCE

$HC1

$BR1
NULL

$TR1
NULL

2.1.3 New data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain))))

    })

ggs

$SUR

$CCE

$HC1

$BR1

$TR1

2.1.3.1 Clumping molecular clocks

2.1.3.1.1 Old data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
    if (nrow(X) > 0)
    ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain))))

    })

ggs

$SUR

$CCE

$HC1

$BR1
NULL

$TR1
NULL

2.1.3.1.2 New data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
  ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %d trees", nrow(X) / length(unique(X$chain))))

    })

ggs

$SUR

$CCE

$HC1

$BR1

$TR1

2.2 Early fossils

2.2.1 New and old data sets

pnts_lks <- readRDS("./output/topology/tree_distance_early_fossils_pairwise_shared.RDS")

ggs <- lapply(pnts_lks, function(X) {

    ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") + 
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain)))) 

    })

ggs

$SUR

$CCE

$HC1

$BR1

$TR1

2.2.2 New data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))
  })

ggs

$SUR

$CCE

$HC1

$BR1

$TR1

2.2.3 Old data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  
  if (nrow(X) > 0)
  ggplot(data = X, aes(x = x, y = y, fill = generation)) +
    geom_path(alpha = 0.25, aes(colour = generation), size = 0.75, show.legend = FALSE) +
    scale_colour_gradient(low = "red", high = "yellow") +
    geom_point(shape = 21, size = 4, colour = "white", alpha = 0.6) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_fill_gradientn(colours = mako(256)) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))
  })

ggs

$SUR

$CCE

$HC1

$BR1
NULL

$TR1
NULL

2.2.4 Clumping molecular clocks

2.2.4.1 Old data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("old", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
    if (nrow(X) > 0)
    ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))

    })

ggs

$SUR

$CCE

$HC1

$BR1
NULL

$TR1
NULL

2.2.4.2 New data set

ggs <- lapply(pnts_lks, function(X) {
     
  X <- X[grep("new", X$chain), ]   
  X$clock <- ifelse(grepl("global", X$chain), "Global", "Uexp")
  X$chain <- gsub("_global|_Uexp", "", X$chain)
  
  ggplot(data = X, aes(x = x, y = y)) +
    geom_path(alpha = 0.1, size = 0.75, show.legend = FALSE) +
    geom_point(aes(shape = clock,  col = clock), size = 4) +
      theme(panel.background = element_blank(), axis.line = element_line(color = "grey"),
      panel.spacing = unit(0.1, "lines"), axis.title.x = element_text(vjust = -0.5), axis.title.y = element_text(vjust = 1.5)) +
    facet_wrap(~chain, nrow = 4) +
      scale_color_viridis_d(option = "G", begin = 0.2, end = 0.8, alpha = 0.75) +
      labs(x = "MDS_v1", y = "MDS_v2") +
    ggtitle(label = sprintf("Tree space for %1.0f trees", nrow(X) / length(unique(X$chain))))

    })

ggs

$SUR

$CCE

$HC1

$BR1

$TR1

3 Topological space dissimilarity

3.1 Observations per lek and alignment

All fossils

# similarity
out <- lapply(pnts_lks, function(x){
  
  Y <- space_similarity(formula = chain ~ x + y, data = x, pb = FALSE, type = "density.overlap")
  
  Y$mean.overlap <- Y$mean.overlap #-  mean(Y$mean.overlap)
  # Y$overlap <- Y$overlap / max(Y$overlap)
  
  return(Y)  
})

topo_similarity_all <- do.call(rbind, out)

topo_similarity_all2 <- topo_similarity_all[, c("group.1",  "group.2", "overlap.1in2")]
names(topo_similarity_all2) <- c("group.1", "group.2", "mean.overlap")

topo_similarity_all3 <- topo_similarity_all[, c("group.2",  "group.1", "overlap.2in1")]
names(topo_similarity_all3) <- c("group.2", "group.1", "mean.overlap")

topo_similarity_all <- rbind(topo_similarity_all2, topo_similarity_all3)

topo_similarity_all$alignment <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][1])

topo_similarity_all$data.set <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][2])

topo_similarity_all$fossils <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][3])

topo_similarity_all$models <- sapply(topo_similarity_all$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][4])


topo_similarity_all$lek <- substr(rownames(topo_similarity_all), 0 , 3)

tab <- as.data.frame(table(topo_similarity_all$lek, topo_similarity_all$alignment))
names(tab) <- c("lek", "model", "count")

tab

lek	model	count
BR1	all.equal	18
CCE	all.equal	76
HC1	all.equal	76
SUR	all.equal	76
TR1	all.equal	18
BR1	optimal	10
CCE	optimal	44
HC1	optimal	44
SUR	optimal	44
TR1	optimal	10
BR1	prank	2
CCE	prank	12
HC1	prank	12
SUR	prank	12
TR1	prank	2

Early fossils

out <- lapply(pnts_lks, function(x){
  
  Y <- space_similarity(formula = chain ~ x + y, data = x, pb = FALSE, type = "mean.density.overlap")
  
  Y$mean.overlap <- Y$mean.overlap -  mean(Y$mean.overlap)
  # Y$overlap <- Y$overlap / max(Y$overlap)
  
  return(Y)  
})

topo_similarity_early <- do.call(rbind, out)

topo_similarity_early2 <- topo_similarity_early[, c("group.1",  "group.2", "overlap.2in1")]
names(topo_similarity_early2) <- c("group.2", "group.1", "mean.overlap")

topo_similarity_early3 <- topo_similarity_early[, c("group.1",  "group.2", "overlap.1in2")]
names(topo_similarity_early3) <- c("group.2", "group.1", "mean.overlap")

topo_similarity_early <- rbind(topo_similarity_early2, topo_similarity_early3)

topo_similarity_early$alignment <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][1])

topo_similarity_early$data.set <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][2])

topo_similarity_early$fossils <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][3])

topo_similarity_early$models <- sapply(topo_similarity_early$group.1, function(x) strsplit(x, "_", fixed = TRUE)[[1]][4])

topo_similarity_early$lek <- substr(rownames(topo_similarity_early), 0 , 3)

tab <- as.data.frame(table(topo_similarity_early$lek, topo_similarity_early$alignment))
names(tab) <- c("lek", "model", "count")

tab

lek	model	count
BR1	all.equal	2
CCE	all.equal	12
HC1	all.equal	12
SUR	all.equal	12
TR1	all.equal	2
BR1	optimal	10
CCE	optimal	44
HC1	optimal	44
SUR	optimal	44
TR1	optimal	10
BR1	prank	18
CCE	prank	76
HC1	prank	76
SUR	prank	76
TR1	prank	18

3.2 Plot dissimilarity by treatment

topo_similarity_all$fossils <- "all"
topo_similarity_early$fossils <- "early"

shared_columns <- intersect(names(topo_similarity_all), names(topo_similarity_early))

topo_similarity <- rbind(topo_similarity_all[, shared_columns], topo_similarity_early[, shared_columns])

topo_similarity$dissimilarity <- 1 - topo_similarity$mean.overlap

# raincloud plot:
ggplot(topo_similarity,
       aes(x = alignment, y = dissimilarity, fill = alignment)) +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    # fill = fill_color,
    alpha = 0.5,
    # custom bandwidth
    adjust = .5,
    # adjust height
    width = .6,
    .width = 0,
    # move geom to the cright
    justification = -.2,
    point_colour = NA
  ) +
  geom_boxplot(# fill = fill_color,
    width = .15,
    # remove outliers
    outlier.shape = NA) +
    # add justified jitter from the {gghalves} package
    gghalves::geom_half_point(
      # color = fill_color,
      # draw jitter on the left
      side = "l",
      # control range of jitter
      range_scale = .4,
      # add some transparency
      alpha = .5,
      transformation = ggplot2::position_jitter(height = 0)
    ) +
      scale_fill_viridis_d(option = "G",
                           end = 0.8,
                           alpha = 0.6) +
           scale_x_discrete(labels = c(
        "natural" = "Natural",
        "artificial" = "Artificial",
        "simulated" = "Simulated"
      )) +
      theme(legend.position = "none") +
      facet_wrap(~ data.set + models + fossils, ncol = 2) +
  labs(x = "Alignment", y = "Topologic space dissimilarity") 

3.3 Stats

3.3.1 Dissimilarity of the topological space

# topo_similarity$dissimilarity[topo_similarity$dissimilarity <= 0] <- 0.001 

topo_similarity$alignment <- factor(topo_similarity$alignment, levels = c("all.equal", "optimal", "prank"))

topo_similarity$data.set <- factor(topo_similarity$data.set, levels = c("new", "old"))

topo_similarity$models <- factor(topo_similarity$models, levels = c("Uexp", "global"))

topo_similarity$fossils <- factor(topo_similarity$fossils, levels = c("early", "all"))

mod <- brm(
  formula = dissimilarity ~ alignment + data.set + models + fossils + (1 | lek), 
  data = topo_similarity,
  iter = iterations,
  family = gaussian(),
  chains = chains,
  cores = chains,
   control = list(adapt_delta = 0.99, max_treedepth = 15),
  file = "./data/processed/regression_results_topological_similiarity",
  file_refit = "always")

mod <- add_criterion(mod, criterion = c("loo"))

null_mod <- brm(
  formula = dissimilarity ~ 1, 
  data = topo_similarity,
  iter = iterations,
  family = gaussian(),
  chains = chains,
  cores = chains,
   control = list(adapt_delta = 0.99, max_treedepth = 15),
  file = "./data/processed/null_regression_results_topological_similiarity",
  file_refit = "always"
)

null_mod <- add_criterion(null_mod, criterion = c("loo"))

Model selection:

mod <- readRDS("./data/processed/regression_results_topological_similiarity.rds")

null_mod <- readRDS("./data/processed/null_regression_results_topological_similiarity.rds")

loo_diff <- loo::loo_compare(mod, null_mod)

as.data.frame(loo_diff[,1:2], row.names = c("model", "null model"))

	elpd_diff	se_diff
model	0.00000	0.000000
null model	-33.31634	8.978687

extended_summary(fit = mod,
    n.posterior = 1000,
    fill =  mako(10)[7],
    trace.palette = viridis,
    remove.intercepts = TRUE,
    highlight = TRUE,
    trace = TRUE,
    return = FALSE, 
    print.name = FALSE
)

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	Intercept-student_t(3, 0.9, 2.5) sd-student_t(3, 0, 2.5) sigma-student_t(3, 0, 2.5)	dissimilarity ~ alignment + data.set + models + fossils + (1 | lek)	5000	4	1	2500	2 (2e-04%)	0	6259.166	5572.907	2016324370

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
b_alignmentoptimal	0.003	-0.057	0.063	1.001	6969.127	7503.027
b_alignmentprank	0.125	0.054	0.197	1	6259.166	6972.068
b_data.setold	0.188	0.138	0.237	1	8340.606	6680.098
b_modelsglobal	0.011	-0.035	0.056	1	8505.640	5572.907
b_fossilsall	0.090	0.029	0.149	1	6927.156	6884.685

Contrasts for alignment strategies:

# contrasts
contrasts(
    fit = mod,
    sort.levels = c("prank", "optimal", "all.equal"),
    predictor = "alignment",
    n.posterior = 2000,
    level.sep = " VS ",
    html.table = TRUE,
    plot = TRUE,
    highlight = TRUE,
    fill = mako(10)[7]
)

	Contrasts	Estimate	Est.Error	l-95% CI	u-95% CI
1	prank VS optimal	0.122	0.031	0.061	0.183
2	prank VS all.equal	0.125	0.037	0.054	0.197
3	optimal VS all.equal	0.003	0.031	-0.057	0.063

custom_ppc(fit = mod, group = "alignment")

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggridges_0.5.6       cowplot_1.1.3        brmsish_1.0.0       
 [4] brms_2.22.0          Rcpp_1.0.13          PhenotypeSpace_0.1.0
 [7] ape_5.8              ggplot2_3.5.1        kableExtra_1.4.0    
[10] knitr_1.48           viridis_0.6.5        viridisLite_0.4.2   
[13] pbapply_1.7-2       

loaded via a namespace (and not attached):
  [1] tensorA_0.36.2.1       rstudioapi_0.16.0      jsonlite_1.8.9        
  [4] magrittr_2.0.3         TH.data_1.1-2          estimability_1.5.1    
  [7] spatstat.utils_3.1-0   farver_2.1.2           rmarkdown_2.28        
 [10] vctrs_0.6.5            spatstat.explore_3.3-2 terra_1.7-78          
 [13] htmltools_0.5.8.1      curl_5.2.2             distributional_0.5.0  
 [16] tidybayes_3.0.7        raster_3.6-26          StanHeaders_2.32.10   
 [19] KernSmooth_2.23-24     htmlwidgets_1.6.4      plyr_1.8.9            
 [22] sandwich_3.1-0         emmeans_1.10.3         zoo_1.8-12            
 [25] igraph_2.0.3           lifecycle_1.0.4        pkgconfig_2.0.3       
 [28] Matrix_1.7-0           R6_2.5.1               fastmap_1.2.0         
 [31] digest_0.6.37          colorspace_2.1-1       CircStats_0.2-6       
 [34] tensor_1.5             vegan_2.6-8            labeling_0.4.3        
 [37] fansi_1.0.6            spatstat.sparse_3.1-0  adehabitatHR_0.4.22   
 [40] polyclip_1.10-7        abind_1.4-8            mgcv_1.9-1            
 [43] compiler_4.4.1         proxy_0.4-27           withr_3.0.1           
 [46] backports_1.5.0        inline_0.3.19          DBI_1.2.3             
 [49] QuickJSR_1.4.0         pkgbuild_1.4.4         MASS_7.3-61           
 [52] classInt_0.4-10        loo_2.8.0              permute_0.9-7         
 [55] tools_4.4.1            units_0.8-5            goftest_1.2-3         
 [58] quadprog_1.5-8         glue_1.8.0             nlme_3.1-165          
 [61] grid_4.4.1             sf_1.0-17              checkmate_2.3.2       
 [64] reshape2_1.4.4         cluster_2.1.6          ade4_1.7-22           
 [67] generics_0.1.3         gtable_0.3.5           spatstat.data_3.1-2   
 [70] class_7.3-22           tidyr_1.3.1            sp_2.1-4              
 [73] xml2_1.3.6             utf8_1.2.4             spatstat.geom_3.3-2   
 [76] pillar_1.9.0           ggdist_3.3.2           stringr_1.5.1         
 [79] posterior_1.6.0        splines_4.4.1          dplyr_1.1.4           
 [82] lattice_0.22-6         gghalves_0.1.4         survival_3.7-0        
 [85] deldir_2.0-4           tidyselect_1.2.1       arrayhelpers_1.1-0    
 [88] gridExtra_2.3          V8_4.4.2               svglite_2.1.3         
 [91] stats4_4.4.1           xfun_0.48              adehabitatMA_0.3.17   
 [94] bridgesampling_1.1-2   matrixStats_1.4.1      adehabitatLT_0.3.28   
 [97] rstan_2.32.6           stringi_1.8.4          yaml_2.3.10           
[100] boot_1.3-30            evaluate_1.0.0         codetools_0.2-20      
[103] tibble_3.2.1           cli_3.6.3              RcppParallel_5.1.9    
[106] xtable_1.8-4           systemfonts_1.1.0      munsell_0.5.1         
[109] spatstat.random_3.3-1  coda_0.19-4.1          svUnit_1.0.6          
[112] spatstat.univar_3.0-1  parallel_4.4.1         rstantools_2.4.0      
[115] bayesplot_1.11.1       Brobdingnag_1.2-9      phangorn_2.11.1       
[118] mvtnorm_1.3-1          scales_1.3.0           e1071_1.7-16          
[121] purrr_1.0.2            rlang_1.1.4            fastmatch_1.1-4       
[124] multcomp_1.4-25